################################################################################
##############################06.Fixed Effects##################################
################################################################################

#####Load data from Data_Generation_CPS#####

pols <- readRDS('pols_bjps')

#####Load plm package#####

library(plm)
library(pcse)
library(lmtest)
library(dplyr)

#####Table 1: Model 1#####

fixed1 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc1 + 
               scale_socioeconomicsal +
               scale_galtansdw +
               scale_enep +
                scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(fixed1, vcov = vcovHC)

fixed2 <- plm(scale_econsdw ~ scale_dispgini + 
                scale_ggini_inc2 + 
                scale_socioeconomicsal+
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(fixed2, vcov = vcovHC)

fixed3 <- plm(scale_econsdw ~ scale_dispgini + 
                scale_ggini_inc3 + 
                scale_socioeconomicsal +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(fixed3, vcov = vcovHC)

fixed4 <- plm(scale_econsdw ~ scale_dispgini + 
                scale_ggini_inc4 + 
                scale_socioeconomicsal +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(fixed4, vcov = vcovHC)

fixed5 <- plm(scale_econsdw ~ scale_dispgini + 
                scale_ggini_inc5 + 
                scale_socioeconomicsal +
                scale_galtansdw +
                scale_enep +
                scale_unemp,
              data = pols, 
              na.action = na.omit, 
              model = "within")

coeftest(fixed5, vcov = vcovHC)


#####Table 1: Model 1 Output#####

#Model fit statistics
fit1b <- c(summary(fixed1)$r.squared, summary(fixed1)$fstatistic$statistic, summary(fixed1)$fstatistic$p.value)
fit2b <- c(summary(fixed2)$r.squared, summary(fixed2)$fstatistic$statistic, summary(fixed2)$fstatistic$p.value)
fit3b <- c(summary(fixed3)$r.squared, summary(fixed3)$fstatistic$statistic, summary(fixed3)$fstatistic$p.value)
fit4b <- c(summary(fixed4)$r.squared, summary(fixed4)$fstatistic$statistic, summary(fixed4)$fstatistic$p.value)
fit5b <- c(summary(fixed5)$r.squared, summary(fixed5)$fstatistic$statistic, summary(fixed5)$fstatistic$p.value)


fit.statsb <- rbind(fit1b, fit2b, fit3b, fit4b, fit5b)

mean(fit.statsb[,1]) #R-Squared
mean(fit.statsb[,2]) #Adjusted R-Squred
mean(fit.statsb[,3]) #F-Statistic
mean(fit.statsb[,4]) #p-value of F-Statistic


fe1 <- as.data.frame(summary(fixed1, vcov = vcovHC)$coefficients) 

fe2 <- as.data.frame(summary(fixed2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fe3 <- as.data.frame(summary(fixed3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fe4 <- as.data.frame(summary(fixed4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fe5 <- as.data.frame(summary(fixed5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fe <- cbind(fe1, fe2, fe3, fe4, fe5)

dot.data.fe <- dot.data.fe[ , order(names(dot.data.fe))]

dot.table.fe <- dot.data.fe %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fe, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fe, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fe, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fe, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fe$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience","GALTAN Polarization", "ENEP",  "Unemployment")

dot.table.fe$Predictor <- factor(dot.table.fe$Predictor, 
                               levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","ENEP",  "GALTAN Polarization", "Unemployment"))

dot.table.fe$Predictor <- factor(dot.table.fe$Predictor, levels=rev(levels(dot.table.fe$Predictor)))

rownames(dot.table.fe) <- dot.table.fe$Predictor
dot.table.fe <- dot.table.fe[,-7]

print(xtable::xtable(dot.table.fe, digits = 3))



#####Table 1: Model 2#####

gg1 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc1 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc1 + 
             scale_galtansdw +
             scale_enep +
             scale_unemp,
           data = pols,
           na.action = na.omit,
           model = 'within')

coeftest(gg1, vcov = vcovHC)

gg2 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc2 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc2 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg2, vcov = vcovHC)

gg3 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc3 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc3 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg3, vcov = vcovHC)

gg4 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc4 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc4 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg4, vcov = vcovHC)

gg5 <- plm(scale_econsdw ~ scale_dispgini + 
             scale_ggini_inc5 + 
             scale_socioeconomicsal +
             scale_dispgini : scale_ggini_inc5 + 
             scale_galtansdw +
             scale_enep+
             scale_unemp,
           data = pols, 
           na.action = na.omit, 
           model = "within")

coeftest(gg5, vcov = vcovHC)

#####Table 1: Model 2 Output#####

fu1 <- as.data.frame(summary(gg1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(gg2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(gg3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(gg4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(gg5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini", "Gini*Income Dif.", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))


#Model Fit Statistics#

fit1i <- c(summary(gg1)$r.squared, summary(gg1)$fstatistic$statistic, summary(gg1)$fstatistic$p.value)
fit2i <- c(summary(gg2)$r.squared, summary(gg2)$fstatistic$statistic, summary(gg2)$fstatistic$p.value)
fit3i <- c(summary(gg3)$r.squared, summary(gg3)$fstatistic$statistic, summary(gg3)$fstatistic$p.value)
fit4i <- c(summary(gg4)$r.squared, summary(gg4)$fstatistic$statistic, summary(gg4)$fstatistic$p.value)
fit5i <- c(summary(gg5)$r.squared, summary(gg5)$fstatistic$statistic, summary(gg5)$fstatistic$p.value)


fit.statsi <- rbind(fit1i, fit2i, fit3i, fit4i, fit5i)

print(mean(fit.statsi[,1])) #R-Squared
print(mean(fit.statsi[,2])) #Adjusted R-Squred
print(mean(fit.statsi[,3])) #F-Statistic
print(mean(fit.statsi[,4])) #p-value of F-Statistic

#####Table 1: Model 3#####

sal1 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc1 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal1, vcov = vcovHC)

sal2 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc2 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal2, vcov = vcovHC)

sal3 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc3 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal3, vcov = vcovHC)

sal4 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc4 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal4, vcov = vcovHC)

sal5 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc5 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_socioeconomicsal + 
              scale_galtansdw +
              scale_enep +
              scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(sal5, vcov = vcovHC)

#####Table 1: Model 3 Output#####

fu1 <- as.data.frame(summary(sal1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(sal2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(sal3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(sal4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(sal5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Salience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini", "Gini*Salience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))

#Model Fit Statistics#

fit1h <- c(summary(sal1)$r.squared, summary(sal1)$fstatistic$statistic, summary(sal1)$fstatistic$p.value)
fit2h <- c(summary(sal2)$r.squared, summary(sal2)$fstatistic$statistic, summary(sal2)$fstatistic$p.value)
fit3h <- c(summary(sal3)$r.squared, summary(sal3)$fstatistic$statistic, summary(sal3)$fstatistic$p.value)
fit4h <- c(summary(sal4)$r.squared, summary(sal4)$fstatistic$statistic, summary(sal4)$fstatistic$p.value)
fit5h <- c(summary(sal5)$r.squared, summary(sal5)$fstatistic$statistic, summary(sal5)$fstatistic$p.value)


fit.statsh <- rbind(fit1h, fit2h, fit3h, fit4h, fit5h)

print(mean(fit.statsh[,1])) #R-Squared
print(mean(fit.statsh[,2])) #Adjusted R-Squred
print(mean(fit.statsh[,3])) #F-Statistic
print(mean(fit.statsh[,4])) #p-value of F-Statistic


#####Table 1: Model 4#####

full1 <- plm(scale_econsdw ~ scale_dispgini + 
              scale_ggini_inc1 + 
              scale_socioeconomicsal +
              scale_dispgini : scale_ggini_inc1 + 
              scale_dispgini:scale_socioeconomicsal +
              scale_ggini_inc1: scale_socioeconomicsal +
              scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
              scale_galtansdw +
              scale_enep +
               scale_unemp,
            data = pols, 
            na.action = na.omit, 
            model = "within")

coeftest(full1, vcov = function(x) vcovHC(x))

full2 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc2 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc2 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc2: scale_socioeconomicsal +
               scale_ggini_inc2:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full2, vcov = vcovHC)

full3 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc3 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc3 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc3: scale_socioeconomicsal +
               scale_ggini_inc3:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full3, vcov = vcovHC)

full4 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc4 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc4 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc4: scale_socioeconomicsal +
               scale_ggini_inc4:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full4, vcov = vcovHC)

full5 <- plm(scale_econsdw ~ scale_dispgini + 
               scale_ggini_inc5 + 
               scale_socioeconomicsal +
               scale_dispgini : scale_ggini_inc5 + 
               scale_dispgini*scale_socioeconomicsal +
               scale_ggini_inc5: scale_socioeconomicsal +
               scale_ggini_inc5:scale_socioeconomicsal : scale_dispgini +
               scale_galtansdw +
               scale_enep +
               scale_unemp,
             data = pols, 
             na.action = na.omit, 
             model = "within")

coeftest(full5, vcov = vcovHC)

##### Table 1: Model 4 Output#####

fit1c <- c(summary(full1)$r.squared, summary(full1)$fstatistic$statistic, summary(full1)$fstatistic$p.value)
fit2c <- c(summary(full2)$r.squared, summary(full2)$fstatistic$statistic, summary(full2)$fstatistic$p.value)
fit3c <- c(summary(full3)$r.squared, summary(full3)$fstatistic$statistic, summary(full3)$fstatistic$p.value)
fit4c <- c(summary(full4)$r.squared, summary(full4)$fstatistic$statistic, summary(full4)$fstatistic$p.value)
fit5c <- c(summary(full5)$r.squared, summary(full5)$fstatistic$statistic, summary(full5)$fstatistic$p.value)


fit.statsc <- rbind(fit1c, fit2c, fit3c, fit4c, fit5c)

mean(fit.statsc[,1]) #R-Squared
mean(fit.statsc[,2]) #Adjusted R-Squred
mean(fit.statsc[,3]) #F-Statistic
mean(fit.statsc[,4]) #p-value of F-Statistic

fu1 <- as.data.frame(summary(full1, vcov = vcovHC)$coefficients) 

fu2 <- as.data.frame(summary(full2, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate2 = Estimate, 
                Std.Error2 = `Std. Error`, 
                `t-value2` = `t-value`, 
                `Pr(>|t|)2` = `Pr(>|t|)`)
fu3 <- as.data.frame(summary(full3, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate3 = Estimate, 
                Std.Error3 = `Std. Error`, 
                `t-value3` = `t-value`, 
                `Pr(>|t|)3` = `Pr(>|t|)`)
fu4 <- as.data.frame(summary(full4, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate4 = Estimate, 
                Std.Error4 = `Std. Error`, 
                `t-value4` = `t-value`, 
                `Pr(>|t|)4` = `Pr(>|t|)`)
fu5 <- as.data.frame(summary(full5, vcov = vcovHC)$coefficients) %>%
  dplyr::rename(Estimate5 = Estimate, 
                Std.Error5 = `Std. Error`, 
                `t-value5` = `t-value`, 
                `Pr(>|t|)5` = `Pr(>|t|)`)

dot.data.fu <- cbind(fu1, fu2, fu3, fu4, fu5)

dot.data.fu <- dot.data.fu[ , order(names(dot.data.fu))]

dot.table.fu <- dot.data.fu %>%
  transmute(pe = rowMeans(dplyr::select(dot.data.fu, Estimate, Estimate2, Estimate3, Estimate4, Estimate5)),
            se = rowMeans(dplyr::select(dot.data.fu, `Std. Error`, Std.Error2, Std.Error3, Std.Error4, Std.Error4, Std.Error5)),
            t.value = rowMeans(dplyr::select(dot.data.fu, `t-value`, `t-value2`, `t-value3`, `t-value4`, `t-value5`)), 
            p.value = rowMeans(dplyr::select(dot.data.fu, `Pr(>|t|)`, `Pr(>|t|)2`, `Pr(>|t|)3`, `Pr(>|t|)4`, `Pr(>|t|)5`))) %>%
  mutate(lwr = pe - se*1.96) %>%
  mutate(up = pe + se*1.96)

dot.table.fu$Predictor <- c("Disposable Gini", "Income Dif.", "Economic Salience", "GALTAN Polarization", "ENEP", "Unemployment", "Gini*Income Dif.", "Gini*Salience", "Income Dif.*Salience", "Gini*Income Dif.*Saience")

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, 
                                 levels = c( "Income Dif.", "Economic Salience", "Disposable Gini","Income Dif.*Salience", "Gini*Income Dif.", "Gini*Salience","Gini*Income Dif.*Saience", "GALTAN Polarization","ENEP", "Unemployment"))

dot.table.fu$Predictor <- factor(dot.table.fu$Predictor, levels=rev(levels(dot.table.fu$Predictor)))

library(ggplot2)

#dot.table.fu %>% ggplot(aes(x = Predictor, y = pe)) + 
#  geom_point() + 
#  geom_hline(aes(yintercept = 0), linetype = 'dashed') +
#  geom_linerange(aes(ymin = lwr, ymax = up))+
#  coord_flip() +
#  labs(title = "Point Estimates of Model Regressing Predictors #on Economic Polarization (n = 83)", y = "Point Estimates (95% #Confidence Interval)", x = "Predictor")+
#  theme_bw()

rownames(dot.table.fu) <- dot.table.fu$Predictor
dot.table.fu <- dot.table.fu[,-7]

print(xtable::xtable(dot.table.fu, digits = 3))


#####Do I need the fixed effects?#####


form <- scale_econsdw ~ scale_dispgini + 
  scale_ggini_inc1 + 
  scale_socioeconomicsal +
  scale_dispgini : scale_ggini_inc1 + 
  scale_dispgini*scale_socioeconomicsal +
  scale_ggini_inc1: scale_socioeconomicsal +
  scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
  scale_galtansdw +
  scale_enep +
  scale_unemp

#robust hausman used to adjust the v-cov matrix appropriately
phtest(form, data = pols, method = "aux", vcov = vcovHC) #Definitely need them#


#####Do I need fixed time effects too?#####

fixedtime<- plm(scale_econsdw ~ scale_dispgini + 
                  scale_ggini_inc1 + 
                  scale_socioeconomicsal +
                  scale_dispgini : scale_ggini_inc1 + 
                  scale_dispgini*scale_socioeconomicsal +
                  scale_ggini_inc1: scale_socioeconomicsal +
                  scale_ggini_inc1:scale_socioeconomicsal : scale_dispgini +
                  scale_galtansdw +
                  scale_enep +
                  scale_unemp +
                  factor(electionyear),
              data = pols, 
              na.action = na.omit,
              model = 'within')

pFtest(fixedtime, full1)#No evidence for two-way effects#

#####Do I need to worry about Serial Correlation#####

plmtest(full1, c("time"), type=("bp"))#No evidence of time effects#

pwartest(full1) #No significant evidence of serial correlation#


